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An extension of the the density matrix renormalization group 
(DMRG) method is presented. Besides the two groups or classes of 
block states considered in White's formulation, the retained m states 
and the neglected ones, we introduce an intermediate group of block 
states having the following p largest eigenvalues Ai of the reduced den- 
sity matrix: Ai > . . . A m > A m +i > . . . > A m + P . These states are taken 
into account when they contribute to intrablock transitions but are ne- 
glected when they participate in more delocalized interblock fluctua- 
tions. Applications to one-dimensional models (Heisenberg, Hubbard 
and dimerized tight-binding) show that in this way the involved com- 
puter resources can be reduced without significant loss of accuracy. The 
efficiency and accuracy of the method is analyzed by varying m and p 
and by comparison with standard DMRG calculations. A Hamiltonian- 
independent scheme for choosing m and p and for extrapolating to the 
limit where m and p are infinite is provided. Finally, an extension of 
the 3-classes approach is outlined, which incorporates the fluctuations 
between the p states of different blocks as a self-consistent dressing of 
the block interactions among the retained m states. 



I. INTRODUCTION 

The density matrix renormalization group (DMRG) method introduced by 
S. R. White in 1992EI has rapidly become one of the leading numerical tools for 
the study of one-dimensional (ID) and quasi-lD strongly correlated systems. Ac- 
curacy and flexibility are two main qualities of this approach. Indeed, the DMRG 
method allows to study large many-body systems containing hundreds of sites and 
to determine several ground-state and low-lying excited-state properties with a pre- 
cision nearly comparable to that of conventional exact diagonalization methods. 
The success of this technique is readily demonstrated Jpy its wide range of applica- 
tions. These concern Heisenberg spin HamiltoniansaB, the t-J modeH, Hubbard- 
like Hamiltonians on various lattices (e.g., dimerized polymjsr chains and Bcthc 
lattices)!] 0, the Kondo impurity problemQ, acoustic phononsEI, etc. 

In its simplest form, the method can be regarded as an iterative projection tech- 
nique, which allows to include the most relevant part of the ground-state wave- 
function on a limited number of many-body states and which leads to renormalized 
block interactions between different regions in real space. The accuracy of the cal- 
culation is mainly controlled by the number m of basis states which are retained in 
each block after a renormalization iteration. When combined with efficient diago- 
nalization techniques for sparse matrices (e.g., Lanczos or Davidson methods) the 
ground-state energy of an infinite ID system can be extrapolated with a remarkable 
precision at a much lower computational cost than conventional exact diagonaliza- 
tions (e.g., straightforward Lanczos). However, as the systems under study become 
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more complex (e.g., itinerant electrons in rings, ladders or multi-band wires) and 
as more delicate -properties are considered (e.g, charge and spin gaps, dynamical re- 
sponse functional, etc.) the DMRG method faces its own limitations, particularly 
since accuracy improves slowly with increasing m. In such cases the value of m and 
thus the size of the projected Hilbert space required for an appropriate accuracy 
become excessively large. The involved computer resources exceed reasonable lim- 
its and in many relevant situations the calculations are simply non-feasible. It is 
therefore of considerable interest to improve the efficiency of the DMRG method in 
order to reduce the computational effort without significant loss of accuracy. 

In this paper an extension of DMRG method is presented which takes better 
profit from the hierarchy among the states spanning the renormalized blocks. As 
proposed by White in his original papeiu each block b is best described in terms 
of the eigenstates \i)b of the reduced ground-state density matrix. The correspond- 
ing eigenvalues A; represent the population of the block-state \i)b and thus pro- 
vide a natural hierarchy among them. Two classes of states are considered in 
the original DMRG procedure: the m states having the largest Xi which are re- 
tained and the remaining states which are completely neglected. In this paper 
a larger flexibility is introduced by distinguishing an additional group of p states 
(Ai > . . . > A m > A m+ i > . . . Xm+p) which are primarily important for the de- 
scription of the ground-state wave function at block b and which can be neglected 
when they participate to long-range interblock fluctuations. It is shown that in this 
way the computer resources involved in the calculations can be appreciably reduced 
without introducing significant additional inaccuracies. Moreover, it is also inter- 
esting to vary the values of m and p and to quantify how this reflects in the physical 
properties in order to gain further insight into the internal many-body structure of 
the renormalized superblocks. 

The rest of the paper is organized as follows. In the next section the 3-classes 
renormalization method is presented and compared with White's original approach!! 
In Sec. [II the method is applied to a few relevant models for which exact or good nu- 
merical results are available. Improvements and limitations of the 3-classes method 
are discussed in particular by analyzing the results as a function of the sizes m 
and p of class-1 and class-2 subspaces. A model independent scheme for optimizing 
the choice of m and p and for extrapolating to the limit where m and p are infi- 
nite is also considered. Finally, Sec. IV summarizes our conclusions and outlines a 
self-consistent extension of the present approach. 



II. 3 CLASSES DMRG METHOD 



In this section we present a density-matrix renormalization method using 3 classes 
of block states. Fiist of all, White's original DMRG algorithm for infinite one- 
dimensional systemsa is briefly recalled for the sake of clarity and in order to facil- 
itate comparisons. 

The DMRG method is an iterative procedure. In the case of infinite open chains 
the ground-state properties of the infinite system are obtained from the extrapola- 
tion of a succession of calculations on finite-length chains. At each renormalization 
group (RG) iteration N, the number of sites L is increased by 2 (L = 2N + 2). 
In spite of the fact that L increases with N, the dimension of the Hilbert space is 
kept constant by means of the following approximation. At each RG iteration the 
Hilbert space of the L-site chain is projected onto a subspace £ spanned by a lim- 
ited number of many-body states ^>ijki which are intended to take into account the 
most relevant part of the ground-state wave function. The ^iju are constructed as 
the antisymmetrized direct product of many-body states, which are usually referred 
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to as "block states." These blocks correspond to different regions in real space as 
illustrated in Fig. [|. Let Tpf (ipf') stand for ith many-body state of the left (right) 
block (1 < i < m) and let cj>j be the jth state of a central (unrenormalized) site 
(1 < j < mi). Then 

*«w = i>i ® 4>j ® 4>k ® i>? ■ (2.1) 

White showed that the error introduced by the truncation of the Hilbert space 
is minimized when the states spanning the subspace of each renormalized block b 
are the eigenvectors corresponding to the m largest eigenvalues of the ground-state 
density matrix reduced to b. In other words, only the block states having the largest 
occupations are kept. 

The DMRG algorithm can be summarized by four main steps involved in going 
from iteration N to iteration N + 1: i) The ground state of the L-site system 
(L = 2N + 2), which corresponds to the Hilbert space S(N) of the iteration N, is 
determined. This involves a calculation in a space of limited dimension (rami) 2 . 
Sparse- matrix diagonalization methods (e.g., Lanczos or Davidson algorithm) are 
normally used, ii) The reduced ground-state density matrices p L and p R of the 
superblocks are calculated. As illustrated in Fig. [|, the superblocks are formed by 
one of the side blocks (L or R) and its neighboring site. The eigenvalues Xf (Af ) 
and eigenvectors ipf (V'f 2 ) oi p L (p ) are determined. Hi) The superblock subspace 
of dimension mmi is projected onto the m most populated states derived in the 
previous step (Ai > . . . > A m ). These are the left and right block-states of the 
next iteration. In this way the central sites have been projected or renormalized 
into the side blocks, iv) The Hilbert space £(N + 1) of the iteration N + 1, whic h 



corresponds to a system with 27V + 4 sites, is constructed according to Eq. ( |2.l|) . 
The renormalized Hamiltonian is determined by performing the appropriate unitary 
transformation and projection. Finally, the procedure loops back to step i) until 
the physical properties extrapolated to L = oo converge. 

Two groups or classes of block states are distinguished in White's DMRG pro- 
cedure. The first one is given by the m most populated states which are retained, 
and the second one by the remaining mm\ — m least populated states which are 
projected out. This criterion certainly respects the hierarchy of states given by the 
eigenvalues or populations of the different block states. However, it results in a 
rather rigid procedure, since the contributions of a block state must be either kept 
in full or completely neglected. Except in cases where the number of states that 
can be kept is large enough to ensure the desired accuracy or when a systematic 
extrapolation to m = oo can be performed, the choice of m is a delicate and drastic 
one. In general there is no clear gap in the populations of the block states which 
could justify where the truncation should take place. It is our purpose to render 
the method more flexible and efficient by the introduction of a new intermediate 
class of block states which are intended to contribute in a restricted way to the 
construction of the Hilbert space of the total system. Physically, one expects that 
the description of the internal degrees of freedom of the superblocks and the renor- 
malization of block-site interactions should be more demanding and more relevant 
to the ground-state properties than the description of fluctuations between left and 
right block-states which are farther apart. The intermediate class of states should 
take this into account. Moreover, with three classes of block states, the transi- 
tion between kept and discarded many-body processes should be smoother and the 
actual choice of m less critical. 

We consider 3 groups or classes of eigenstates of the block density matrix: i) The 
first group (class 1) is given by the m most populated states tpi where 1 < i < m and 
Ai > . . . > A m . ii) The second group (class 2) contains the eigenstates tpi associated 
to the following^ eigenvalues, i.e., m+1 < i < m+p with A m > A„ l+ i > . . . > X m + P - 
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Hi) Finally, the third group (class 3) refers to the remaining mmi — (m +p) states. 
The states of the first and third classes are treated just as in White's algorithm, i.e., 
all direct-product states or configurations of the complete system ^iju which involve 
only class-1 block-states (i, I < m) are taken into account while configurations 
involving a class-3 state (i or I > m + p) are all projected out [see Eq. (2.1)]. 
An intermediate criterion is used for the second class, namely, class-2 states of 
the left superblock are included only in combination with class-1 states of the right 
superblock and vice versa. In other words, direct-product wave- functions f y -M with 
m+l<i<m+p and m+l<l<m+p are neglected. This is illustrated in 

In analogy with configurations interaction methods, the subspace spanned by the 
set of ^ijki with i,l < m can be regarded as a complete active space (CAS), the 
subspace spanned by the ^ijki with i < m and m-\-\<l<m+p or I < m and 
m+l<i<m + p defines then a kind of mono-excited space, and the discarded 
states obtained from with m + 1 < i, I < m + p can be viewed as double 

many-body excitations. Class-2 states have a relatively small weight in the ground- 
state wave-function. One therefore expects that with an appropriate choice of m 
and p the contributions of ^ijki obtained by the tensor product of two such block 
states should be particularly small, eventually even less important than the states 
involving products of a class-1 and class-3 states which are already discarded in 
White's DMRG method. Therefore, the accuracy of a 3-classes renormalization 
should not differ significantly from that of a standard 2-class DMRG-calculations 
retaining m + p states. 

Depending on the values of m and p, the 3-classes procedure may allow an ap- 
preciable reduction of the dimension of the many-body Hilbert space, namely, from 
Dq = m\{m + p) 2 to D = m\(m 2 + 2mp). This is of considerable interest in many 
practical applications which are often limited by computer-memory needs and most 
critically by the computer time required for the determination of the ground-state 
wave-function at each RG iteration. The 3-classes approach can also be applied 
straightforwardly to more complex latticgs such as Bethe lattices, which involve 
renormalizations at more than two blocksQ. In these cases the reduction of compu- 
tational effort is even more important, since Dq/D = (m + p) nh / (m nb +ni,pm nb ~ 1 ) 
where nt, is the number of renormalized blocks (e.g., nj, — z for a Bethe lattice with 
z nearest neighbors). 

The eigenvalues of the density matrix reduced to block b represent the proba- 
bility of block b being in state ty\ (b — L or R). Therefore, the traces Si = YaLi 
S 2 = YTJrn+i \ and s 3 = J2T=,n+ P +i measure the weight of each of the 3 
classes of block states on the ground-state wave function (Si + S 2 + S 3 = 1). Ap- 
proximating for simplicity the norm of the retained part of the system wave-function 
by the product of the traces Sj on the different blocks and keeping only the lead- 
ing terms in S 2 and S3 (S 2 ,S3 <C 1), one estimates the truncation error e(m,p) 
involved in a 3-classes renormalization as e(m,p) — + (™ b )S|, where rib refers 
to the number of renormalized blocks (n& = 2 in the present calculations). If p = 0, 
S 2 = and we recover White's estimation of the truncation error 1 — P m whic h is 
proportional to S3EJ (S3 = 1 — Si = 1 — P m if S 2 = 0). As discussed in Sec. Ill, 
comparison with exact results show that e(m,p) gives very good estimation of the 
errors involved in 3-classes calculations. In practice, e(m,p) can be used either to 
control the accuracy of the results, which often depend on the model Hamiltonian 
and its parameters, or to extrapolate to the limit e(m,p) — > (m,p — > 00) in order 
to further reduce,the errors. A similar procedure has been already applied in the 
standard DMRGa. Moreover, e(m,p) provides a practical criterion for optimizing 
the value of m and p, for instance, when the size of the Hilbert space D(m,p) is 
fixed. In general it is a good choice to keep S| and 2S3 of the same order of mag- 
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nitudc. In fact, for a fixed D(m 1 p) a reduction of £f (e.g, by decreasing p) implies 
an increase of £3 and vice versa. 

III. RESULTS 

The density-matrix renormalization using 3-classes of states has been applied to 
several benchmark problems in order to check the accuracy of the method and to 
analyze the quality of the results as a function of the subspace dimensions.-?/! and 
p. In particular comparison is made with the standard DMRG algorithmic. The 
considered systems include the spin- 1/2 Heisenberg chain, the periodic Hubbard 
chain as a function of correlation strength U /t, and the dimerization of polyacetylene 
using a distance-dependent tight-binding Hamiltonian. In the following our main 
results are discussed. 



A. Spin- 1/2 Heisenberg chain 

In this section we consider the infinite spin-1/2 Heisenberg chain as a repre- 
sentative application of the 3-classes DMRG method to spin systems. This prob- 
lem is particularly interesting from a methodological point of view since the exact 
Bethe-Apjsatz solution as well as extremely good standard DMRG calculations are 
available^. In Table | our results for the ground-state energy E are compared with 
the exact value E ex = — Jin 2 + J/4 and with DMRG calculations using White's 
method which corresponds to p = 0. In all cases, N = 100 renormalization itera- 
tions are performed. This ensures that the extrapolated ground-state energy of the 
infinite-length chain is converged at least within 10 -6 J. Small values of to and p 
are considered in Table | in order to compare the convergence and performance of 
the different methods. Notice, however, that more accurate calculations involving 
a larger number of states (e.g., m = m + p = 44EJ) are not at all demanding with 
present computer facilities. For to and p larger than the values shown in Table |, the 

2- classes and 3-classes DMRG results converge very well to the exact ground-state 
energy and are almost indistinguishable from each other. 

The 3-classes results recover accurately the exact ground-state energy even if as 
few as 5 states are kept in the class-1 subspace or complete active space (CAS). For 
to — 5 and m+p = 21, E — E ex ~ 10~ 4 J, which seems an acceptable accuracy taken 
into account that the size of the Hilbert space is reduced significantly (Dq/D = 2.4). 
Our results show that the intrablock fluctuations between the m class-1 states and 
the p class-2 states play a major role in the description of the ground-state wave 
function. In fact, neglecting them yields so poor results that for m = m+p = 5 the 
DMRG procedure diverges after a few iterations. In contrast, the effects of quantum 
fluctuations between class-2 states of different blocks are much less important for 
the determination for the ground-state energy. These contributions are neglected 
in a 3-classes DMRG calculation with m = 5 (to = 7) and m + p = 21 which 
yields E — E ex = 1.07 x 10~ 4 J (E - E ex = 4.47 x 10~ 5 J) and are included in the 
calculation with m = m + p = 21 which yields E — E ex ~ 2.55 x 10 -5 J. Similar 
conclusions are derived for other values of to and p. 

The calculations for different to and fixed m+p provide further insight on the 
convergence properties of the method. For very small to the difference between 

3- classes and 2-classes results decreases very rapidly for increasing m (see Table | 
for 771 — h and 7). This indicates that a minimum size of the class-1 subspace is 
indispensable for obtaining good accuracy. Further increase of m yields almost the 
same accuracy as the calculation with m = m+p = 21. Nevertheless, the reduction 



5 



of the dimension of the Hilbert space is still appreciable (e.g., Dg/D = 1.4 for 
m = 10). Let us recall that the error ia E decreases slowly with increasing D 
already in the standard DMRG methodEJ. As in the case of very small m, the 
intrablock transitions involving class-2 states are very important for improving the 
accuracy (E — E ex = 2.06 x 10 _4 J for m = m + p = 10). We conclude that the 
3-classes DMRG procedure improves efficiency without significant loss of accuracy, 
provided that the size m of the CAS is not too small (m > 7 in the present case). 

The truncation error e(m,p) = S| + 2S3 provides a very good estimate of the 
accuracy of the calculations for all considered values of m and p. In fact, E — 
E ex follows very closely a linear relation with e(m,p). It is therefore simple to 
extrapolate to the limit e(m,p) — > in order to further improve the results (see 
Table [[J). A similar procedure is used in standard DMRGtl. Moreover, the results 
suggest that the dimensions of class- 1 and class-2 subspaces could be chosen in order 
to optimize accuracy by minimizing e(m,p) for a given computational cost. 

B. Hubbard chain 

The DMRG method converges very rapidly and accurately when applied to prob- 
lems involving strong electron correlations such as the Heisenberg model. However, 
the performance of the method is very sensitive to the strength of the Coulomb 
interactions and to the degree of electron derealization. For instance, in the case 
of the Hubbard chain at half-band filling and for a given number of block states m, 
the relative error Ae in the ground-state energy increases by an order of magnitude 
as we move from the strongly correlated to the uncorrelated limit. For example for 
m = 80 (p = 0), Ae = 8.1 x 10" 5 for U/t = 76, Ae = 4.9 x 10~ 5 for U/t = 4, 
Ae = 3.1 x 10~ 4 for U/t = 1 and Ae = 3.3 x 10" 4 for U/t = 0. Moderate and 
weak electron interactions are thus much more difficult to describe with the DMRG 
method. In addition, for small U/t the algorithm requires a quite larger number 
N of RG steps to converge. Typically, N is about 2-3 times larger for U = than 
for very large U/rB. It is therefore considerably interesting to analyze the accu- 
racy and convergence properties of the 3-classes DMRG method in applications to 
itinerant-electron models such as the Hubbard model. 

In Fig. H results are given for the ground-state energy of the half-filled one- 
dimensional Hubbard chain as function of the Coulomb repulsion strength U/t. 
These were obtained by using the 3-classes DMRG method with various represen- 
tative values of m and p. The standard DMRG algorithm corresponds to p = 0. In 
all cases, the number of renormalization iterations N is increased until the infinite- 
length extrapolated energy has converged within an error smaller than 10~ 6 t. Start- 
ing from the uncorrelated limit, the accuracy improves considerably with increasing 
U/t. Very precise results are obtained particularly for U/t > 4. However, for very 
large U/t (U/t > 16) the relative error may increase slightly. For example, for 
U/t = 76 (U/t = 16) the relative errors Ae = (E - E ex )/E ex are Ae = 4.2 x 10~ 4 
(Ae = 1.1 x 10~ 4 ) for m = m + p = 35, Ae = 8.0 x 10~ 5 (Ae = 3.5 x 10~ 5 ) for 
m = 35 and m+p = 80, and Ae = 8.0 x 10" 5 (Ae = 3.0 x 10" 5 ) for m = m+p = 80. 
Notice that, although the differences between the calculation with m = m + p = 35 
and with m = m +p = 80 are significant, the 3-classes results follow very closely 
the m = m +p = 80 curve. This implies that the largest part of the improvements 
made by increasing m from 35 to 80 are recovered by taking m = 35 and m+p = 80, 
i.e., by including the intrablock fluctuations. Transitions involving class-2 states on 
two different superblocks are not very important in this case. This allows to reduce 
the dimension of the Hilbert space by about 30%. The results seem remarkable, 
particularly if we recall that the 3-classcs method is a controlled approximation in 



6 



the sense that it respects the variational principle. At each RG iteration the energy 
obtained with m = 35 and m + p = 80 is an upper bound for the energy of a 
standard DMRG calculation with m = 80 and p = 0. 

The same trends hold for weaker interactions. For U < At, one observes that the 
differences between the relative errors in the 3-classes and in the standard DMRG 
calculations increase. At the same time, however, the m = m + p = 80 curve 
starts to deviate more significantly from the exact result showing the that block 
states i with smaller eigenvalues of the density matrix (i.e., weaker occupations 
Xi < Xso) should be taken into account. Quantitatively, neglecting the intrablock 
fluctuations involving states with A^ < Aso is a much more important source of error 
than the interblock fluctuations neglected in the 3-classes method [see Fig. §(a)]. 
The additional error introduced by the use of the 3-classes procedure (m = 35 and 
m + p = 80) is only 9% of the total error (m = m + p = 80). Therefore, in order 
to further improve accuracy for U/t < 4 it should be more efficient to increase the 
size p of the class-2 subspace instead of increasing m with p — as in the standard 
DMRG method. This is consistent with our estimation of truncation error which is 
quadratic in the sum of the class-2 eigenvalues £2 and linear in £3. 

In order to investigate this question we compare in Fig. |3|(b) the accuracy of 
the 3-classes approach as a function of the size of class- 1 and class-2 subspaces. 
Several values of m and p are considered which all correspond to approximately 
the same dimension of the Hilbert space of the complete system: D(m, m + p) — 
16[(m +p) 2 -p 2 } ~ L>(80,80) = 102400. Therefore, all these calculations demand 
essentially the same computer resources. One observes that there is an optimum 
choice for m and p for which the relative error (E — E ex )/ E ex is reduced by about 
20-30% with respect to the standard calculation with m = m + p = 80. As in 
the Heisenberg model, a minimum number of class- 1 states is indispensable for 
describing the ground-state of the entire system with a reasonable accuracy. If m is 
too small, for instance m = 15, the results are poorer even if p is very large. Once 
the most important configurations are taken into account as class-1 states, it is more 
useful to maximize p in order to increase the flexibility of the wave-function within 
each superblock rather than to increase the number of class-1 states. This improves 
the treatment of intrablock correlations and allows a more precise renormalization 
of the interactions as the number of sites grows. The results shown in Fig. ||(b) for 
m = 15 {p = 221), m = 20 (p = 170) and m = 35 (p = 109) illustrate this behavior. 

Taking into account the strong dependence of E — E ex on (m,p) and on U/t 
(even for p = 0) it is of considerable practical interest to have a model-independent 
criterion for choosing m and p and for controlling the accuracy of 3-classes DMRG 
calculations. Comparison between Figs. 3(b) and 3(c) shows that the estimation of 
the truncation error e(m,p) = £3 + 2X3 follows the trends in E — E ex . Therefore, 
the minimization of e(m,p) provides an appropriate means of choosing m and p 
for a given dimension of the Hilbert space. Alternatively, one may use m and p 
such that e{m,p) is constant in order to ensure approximately the same accuracy 
in different calculations (e.g., as a function of U/t or of the dimerization 8), Notice, 
that e(m,p) can be easily computed at each RG iteration. 

It is should be also noted that the minimum m, which yields better accuracy than 
the standard DMRG method, as well as the optimal choice for m and p for a given D 
depend on the Coulomb interaction strength U/t. For example, m = 20 and p = 170 
yields the best results for U/t < 1, but for U/t> 2 this is somewhat less accurate 
than the standard DMRG calculation. In the range < U/t < 4 the optimum m 
(p) increases (decreases) as U/t increases. This is the result of two contributions to 
the minimization of Ae, which compete when the dimension of the Hilbert space is 
fixed. Increasing p improves the description of the wave function within each block 
as well as the superblock renormalizations since it reduces £3. However, this is done 
at the expense of decreasing m (or increasing £2) and therefore tends to worsen the 
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description of the interactions between different blocks. For U/t > 4 the results for 
in = 35-80 (p = 1 09-80) are very similar. Close to the optimum (m,p) the energy 
differences remain small. Let us finally recall that the overall convergence and 
accuracy of DMRG methods improve as local charge fluctuations are reduced. The 
quantitative differences between different renormalization or projection strategies 
become therefore less important at large U/t. 

C. Dimerized chains 

The purpose of this section is to discuss the application of the 3-classes DMRG 
method to dimerized polymer chains. The dimerization of polyacetylene is partic- 
ularly interesting since it depends sensitively on the details of the .-wave function 
and thus on the number of states m used in DMRG calculationsH. A distance 
dependent tight-binding Hamiltonian is considered for the ir valence electrons in 
polyacetylene. Repulsive interactions between a electrons are modeled by a pair- 
wise potential E a (rij)B. The hopping integrals t(rij) and the repulsive potential 
E a (rij) are obtained from ab-initio calculations on the gropad-state energy and 
first singlet-triplet excitation energy of the ethylene moleculeta. The main interest 
of the uncorrelated case is methodological, since this is the most difficult limit to 
study using the DMRG method. DMRG calculations for realistic finite values of 
the Coulomb repulsion U are reported in Ref.Q. 

In Fig. |4| results are given for the ground-state energy E(S) of an infinite poly- 
acetylene chain vs. the dimerization S = f^i+i — ^i— i,i |/2, as obtained using the 
model Hamiltonian described above. The exact tight-binding result is also given 
for the sake of comparison. As in the undimerized Hubbard chain, the 3-classes 
calculations with m = 35 and m + p = 150 remove the largest part of the discrep- 
ancies between the m — m + p = 35 and m = m + p = 150 standard calculations. 
In particular the position of the minimum S m i n is improved considerably. However, 
notice that the error AE = E — E ex decreases with increasing 8 and therefore, 
Smin remains overestimated. For m = m + p = 35 (m = m + p = 150) one ob- 
tains S m in — 0.389ao (5 m in — 0.306ao). Using m = 35 and m + p = 150 one finds 
Smin — 0.310ao, while the exact result is S m i n — 0.303ao- As expected, for a given m 
and p, the 3-classes method converges better towards the more demanding 2-classes 
calculation (m + p states per block) as the dimerization increases. 

IV. SUMMARY AND OUTLOOK 

In this paper, an extension of the the density matrix renormalization group 
(DMRG) method has been presented, which achieves larger flexibility and improved 
performances without increasing the computational effort. An intermediate class 
of block states is introduced between the states which are fully taken account and 
those which are discarded at each renormalization step. These states contribute 
only as far as mtrablock fluctuations are concerned and provide a smoother criterion 
for discerning between retained and discarded many-body states. Applications of 
the method to one-dimensional models (Heisenberg, Hubbard and dimerized tight- 
binding) show that in this way the involved computer resources can be reduced 
without significant loss of accuracy. Varying the number p of intermediate states 
relative to the number m of block states which are fully retained one obtains further 
insight into the DMRG method, particularly concerning the internal structure of 
the renormalized blocks and the relative importance of intrablock and interblock 
fluctuations. For example, for a given dimension D(m,p) of the Hilbert space of 
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the renormalized system, optimum values of m and p can be determined which 
minimize the error in the ground-state energy. In these cases the discrepancies with 
exact results are 20-30% smaller than in standard DMRG calculations. Using a 
simple estimation of the truncation error involved in a 3-classes renormalization, a 
model-independent strategy is derived for choosing m and p and for extrapolating 
to the limit where m and p are infinite. In this way reliability and accuracy are 
further improved. 

In order to pursue the development of the 3-classes approach, the fluctuations 
involving class-2 states of different blocks should be taken into account without 
increasing the dimension of the renormalized Hilbert space. This can be done by 
using a self-consistent intermediate Hamiltonian method based on a coupled-cluster 
approximation. Similar techniques are very successful in larae.-scale ab initio config- 
urations interaction (CI) calculations on molecular systemstJ. The model space S 
of the effective Hamiltonian coincides with the Hilbert space of the 3-classes DMRG 
approach (i.e., multiple class-2 states on different blocks are excluded). From the 
calculation of the ground-state wave function "J) within S one obtains relations 
between the wave-function coefficients corresponding to class-1 and class-2 states. 
Assuming a coupled-cluster expansion for |^), these relations allow to infer ap- 
proximately the coefficients of the wave function beyond 5, i.e., when class-2 states 
are present on different blocks. These fluctuation effects can be then included as 
a dressing of the block interactions within S. Since the dressed or intermediate 
Hamiltonian depends on a set of self-consistent equations must be solved at 
each renormalization iteration. Notice, however, that this could be done within the 
Lanczos or Davidson iterative procedure in order to avoid a too large increase of 
computer time. The detailed formulation of this self-consistent approach is rather 
lengthy and beyond the scope of this paper. A summary may be found in the 
Appendix. While an efficient implementation of such an extension should further 
improve performance, it is also true that its complexity risks to limit possible ap- 
plications to very difficult cases where a sufficiently large number of states m and 
p is not practicable or where a systematic extrapolation to the limit where m and 
p are infinite is not possible. 
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APPENDIX 

The aim of this section is to outline a method for taking into account the effect 
of the neglected fluctuations between the class-2 block-states, 

Vi*jki* = ® fa <8> <fe <8> $J = \i*jkl*) (4.1) 

(m < i* ,1* < m + p) without increasing the dimension_pf the Hilbert space. In 
analogy with the theory of intermediate HamiltoniansE^I, we intend to map the 
Hilbert space £ of the standard DMRG procedure (m+p states per block) onto the 
3-classes DMRG Hilbert space S (with m class-1 states and p class-2 states). An 
effective Hamiltonian is thereby defined acting on S as model space, which takes 
into account the effects of the outer space Q = £ — S. 
Let us consider the secular equation of a state ^ a bcd £ S 
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E C m H m , abcd + C abcd E (4.2) 

\ijkl)eS |a/3 7 5)eC 

One would like to rewrite this equation in the following effective Hamiltonian form 

^ t CijklHijkl,abcd + Cabcd^abcd.abcd = C abcd E (4.3) 
\ijkl)eS 

since in this way the contribution of the states belonging to Q is taken into ac- 
count by a diagonal dressing A of the Hamiltonian. A should therefore satisfy the 
equations 

Aabcd,abcd = E -^^—H a p^s,abcd (4-4) 

Notice however that the coefficients C a (3j$ cannot be obtained from the diago- 
nalization of the effective Hamiltonian H e ff = H + A since the latter acts only 
on the model space S. The main difficulty is therefore to find a sound evalua- 
tion of the coefficients Cap^s- This problem can be solved using a coupled-cluster 
approximation. 

Let S L and S R be the cluster operators associated to the left and right renormal- 
ized blocks and let = exp (S L + S R ) be the wave operator going from Sq to £, 
where Sq is generated by the class- 1 blocks states only. S can then be decomposed 
as S = So + Sr + Sl, where Sr (Sl) is the Hilbert space generated out of the 
product of class-2 states in the right (left) renormalized block and class-1 states in 
the left (right) renormalized block. The ground-state ^ G £ can be expanded as 

* = E c w ly w > + ( 4 - 5 ) 
E c ->*w l**i*0 + E c v kl * \^ kl *) + ( 4 - 6 ) 

\i*jkl)£Sn \ijkl*)eSL 

E C i*M* IW> > (4-7) 
\i'jki')eQ 

while the coupled-cluster approximation gives 

* = E c ^ ki + ( 4 - 8 ) 

E E S^ i C ijk i\i*jkl)+ E S^C ijkl \ijkl*) + (4.9) 

E E s^x tl c ijU \i*jkn- (4.io) 



Comparing equations 4J3 and 4J5 one obtains the definition of the cluster operators 

Ci'jkl = E] Si*,iCijkl , 



i<m, Kt 



Cijkl* — E] 



ijkl 



ijkl 
i<m, i<r 



and the evaluation of the unknown coefficients 
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Ci*jM* - ^ S i*,i S ^,i G ijki ■ (4-11) 



These systems of equations are unfortunately over-defined in most practical cases 
[p > to). An optimal choice for S L is obtained by minimizing the square deviation 

C L = ^2 (Ci'jkl - ^ SjjjCjjki) 2 ■ 



i*jkl 
i<i<m+p, !<n 



Thus, 



\ 

-2 CijkiCi'jki — ^ S^i j/CijkiCi'jki 
t ,U J 



where p L refers to the ground-state density matrix reduced to the left renormalized 
block. Finally, one obtains 

S L = {P%fPP£) {Ptp L Pty 1 (4-12) 

where Pf (P^) is the projection operator onto the class-1 (class-2) left block-states. 
Similar equations hold for the right block. 

It should be noted that the clusters operators depend on the ground-state of the 
dressed Hamiltonian, which depends itself on the clusters operators. The whole 
process should then be iterated up to self-consistency. Let us finally recall that self- 
consistent approximations are usually very efficient but do not necessarily satisfy 
the variational principle. 
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TABLE I. Ground-state energy E of the infinite spin-1/2 Heisenberg chain as obtained 
using the 3-classes DMRG method with different numbers of block states m and p in class- 1 
and class-2 subspaces (see Fig. []). E ex = — J In 2 + J/4 refers to the Bethe-ansatz exact 
result and Eq to a standard DMRG calculation with m = 21 block states (p — 0). The 
results labeled m = m+p = oo are obtained by linear extrapolation of the truncation error 
t{m,p) = S| + 2E3 — * using the results for m — 5 and m — 7. Ae = (E — Eo)/(Eo — E ex ) 
and Dq/D — (m +p) 2 /(m 2 + Imp) indicate, respectively, the loss of accuracy and the 
reduction of the dimension of the Hilbert space with respect to standard DMRG (m = 21). 
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FIG. 1. Illustration of the superblock renormalization for infinite one-dimensional 
systems. Crosses represents atomic sites and rectangles renormalized blocks. Notice that 
the total number of sites is increased by 2 at each renormalization iteration N. 
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FIG. 2. Illustration of the different projection subspaces S in White's DMRG methoda 
and in the present 3-classes renormalization. The complete Hilbert space £ involves all 
the antisymmetrized direct products of left and right block-states (see Fig. [l]). The hashed 
areas indicate the subspaces 5 which are kept at each iteration. 
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Fig. 3 
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FIG. 3. DMRG results for the ground-state energy E of the one-dimensional Hubbard 
model at half-band filling as a function of the Coulomb repulsion strength U /t. E ex refers 
to the Bethe-ansatz exact solutiontj. The number m (p) of block states in the class- 1 
(class-2) subspaces is indicated (see Fig. ^). White's DMRG method corresponds to p = 0. 
In (a) the dimensions D(m, m + p) — 16[(m + p) 2 — p 2 ] of the many-body Hilbert spaces 
are D(35, 35) = 19600, £>(35, 80) = 70000, D(80, 80) = 102400. In (b) m and p are varied 
keeping a constant D(m, m + p) ~ D(80, 80). In (c) the estimation of the truncation error 
e(m,p) = E| + 2S3 is shown. 
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Fig. 4 
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FIG. 4. Ground-state energy E of polyacetylene infinite chains (in a.u.) as a function 
of the dimerization S = Iri.j+i — ri_i,j|/2 for R — (r»,i+i + n-i t i)/2 = 2.66oo. The results 
are obtained using a distance-dependent tight-binding model (U-t= 0) which parameters 
are derived from ab initio calculations on the ethylene moleculeO. The number of block 
states m and p kept in the 3-classes DMRG calculations are indicated in the inset. The 
exact tight-binding result is given by the lowest thin curve. 
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